clc,clear
z=readmatrix("E:\roller data\g1-0-1.xlsx"); %读入初始形貌
%z=z*10^3;
[numx,numy]=size(z); %计算初始形貌矩阵的行列数

delta_dg=9.9*10^(-6); %采样间距delta_dg微米（手动输入）

%%计算工件抛光前表面粗糙度
zabs_original=abs(z); %将抛光前的形貌矩阵取绝对值
Ra_origianl=sum(sum(zabs_original))/(numx*numy); %计算抛光前的Ra

Lx=numx*delta_dg; %计算采样区域长度
Ly=numy*delta_dg; %计算采样区域宽度


%%在采用区域内，等间距取numx/numy个采样点
xx=linspace(0,delta_dg*numx*10^6,numx);
yy=linspace(0,delta_dg*numy*10^6,numy);

figure(1)
[xx1,yy1]=meshgrid(xx,yy);  %生成采样点的网格矩阵
surf(xx1,yy1,z,'EdgeColor','interp');  %生成初始表面形貌
zlim([-2 2]); %纵坐标范围[-2,2]
view(-60,80);  %调整视角
title('抛光前的表面形貌');
colorbar;

Vy0=0.06; %输入法向速度，单位m/s；
Vx0=0.12; %输入切向速度，单位m/s
Vx_ave=0.13; %输入平均切向速度，单位m/s
qu=1150*10^6; %输入材料屈服强度，单位Mpa，即：牛顿/平方米
P_ration=[0.3,0.22]; %泊松比，第一项为被加工工件的，第二项为磨粒的
E=[206*10^9,350*10^9]; %弹性模量（被加工工件、磨粒）
R=2*10^(-3); %磨粒半径
A_desnity=3850; %磨粒密度
HV1=680*10^3; %工件的维氏硬度（Gpa）；
KIC=130*10^4; %工件的断裂韧度
m=A_desnity*(4/3)*pi*R^3;  %单个磨粒的质量
V_A=(4/3)*pi*R^3; %计算单个磨粒的体积
Pt=654*10^6*9.8; %定义磨粒在切向方向收到的压力pa，即N/m^2
c3=0.95; %计算单个磨粒横向切割距离时用到的常数，试验测得
n=0.95; %切向方向上磨粒与工件接触面积的指数，用于计算Fx
def=zeros(numx,numy); %定义变形量矩阵
delta_h=zeros(numx,numy); %定义每个峰的去除高度矩阵
V_desnity=0.2624; %输入磨粒的体积浓度（每立方米）
T=180; %抛光时间（分钟）

%%计算法向最大变形深度

%弹性-塑性变形临界速度纵向Vy-ela-pla
Vy_ela_pla=(0.8*0.25^(5/2))^0.5*...
    ((qu^5/A_desnity)*((1-P_ration(1)^2)/E(1)+(1-P_ration(2)^2)/E(2))^4)^0.5*pi*pi;

%弹性-塑性变形临界变形深度
deform_ela_pla=0.25*pi^2*R*qu*((1-P_ration(1)^2)/E(1)+(1-P_ration(2)^2)/E(2))^2;
%弹性-塑性变形临界作用力
F_ela_pla=(2/3)*pi*qu*deform_ela_pla*R;
%塑性-脆性临界变形深度
deform_pla_bri=0.15*(E(1)/HV1)*(KIC/HV1)^2;
%塑性-脆性临界冲击速度
Vy_pla_bri=(Vy_ela_pla^2+qu*(3*deform_pla_bri+deform_ela_pla)*(deform_pla_bri-deform_ela_pla)/(4*R^2*A_desnity))^0.5;

%%计算变形深度（非永久/永久变形）deform/deform_per
if Vy0<=Vy_ela_pla  %弹性阶段

    deform=(0.8*pi*A_desnity*((1-P_ration(1)^2)/E(1)+(1-P_ration(2)^2)/E(2)))^0.4*R*Vy0^0.8;
    deform_per=0; %弹性阶段被加工工件表面不发生永久变形

end

if Vy0>Vy_ela_pla && Vy0<=Vy_pla_bri %塑性阶段
 
    %单个磨粒撞击工件表面产生撞击坑的最大深度
    deform=(1/3)*deform_ela_pla+(1/6)*...
        (16*deform_ela_pla^2+48*R^2*A_desnity*(Vy0^2-Vy_ela_pla^2)*qu^(-1))^0.5;  
    deform_per=deform-deform_ela_pla; %塑性阶段被加工工件表面产生的永久变形深度

end

%%计算单个磨粒切向方向的切割长度
ymax=deform; %宽度方向上被切割截面的最大深度
Dmax=2*sqrt(R^2-(R-deform)^2); %撞击坑最大直径，sqrt(R*deform)为撞击坑的半径
Ax=0.5*ymax*Dmax;
%%c3和n均为实验数据，需要手动输入

%L为单个磨粒的切割长度
L=0.5*c3*Vx0^2/(R^(1.25*n)*ymax^(0.75*n)*Pt); %c1,c2或c3均为实验测得的数据

l_desnity=(3*V_desnity/(4*pi*R^3))^(1/3); %磨粒的线密度
ls=1/l_desnity; %两磨粒的平均间距
Nx=(Lx+L)/ls; %取样区域内某一长度方向上的磨粒数
Ny=(Ly+2*R)/ls; %取样取样内某一宽度方向上的磨粒数
%N=Lx*Ly;
%N=(Lx*Ly)/(ls^2); %单位面积A=Lx*Ly内的磨粒数
N=Nx*Ny;

k = 1;
t_one=k*(ls/Vx_ave)/(60*N); %该表面每次被切割的时间间隔(分钟)
num=T/t_one; %表面被切割的次数

%被第num个磨粒作用后表面形貌的变化
clc,clear
z=readmatrix("E:\roller data\g1-0-1.xlsx"); %读入初始形貌
%z=z*10^3;
[numx,numy]=size(z); %计算初始形貌矩阵的行列数

delta_dg=10*10^(-6); %采样间距delta_dg微米（手动输入）

%%计算工件抛光前表面粗糙度
zabs_original=abs(z); %将抛光前的形貌矩阵取绝对值
Ra_origianl=sum(sum(zabs_original))/(numx*numy); %计算抛光前的Ra

Lx=numx*delta_dg; %计算采样区域长度
Ly=numy*delta_dg; %计算采样区域宽度


%%在采用区域内，等间距取numx/numy个采样点
xx=linspace(0,delta_dg*numx*10^6,numx);
yy=linspace(0,delta_dg*numy*10^6,numy);

figure(1)
[xx1,yy1]=meshgrid(xx,yy);  %生成采样点的网格矩阵
surf(xx1,yy1,z,'EdgeColor','interp');  %生成初始表面形貌
zlim([-2 2]); %纵坐标范围[-2,2]
view(-60,80);  %调整视角
title('抛光前的表面形貌');
colorbar;

Vy0=0.06; %输入法向速度，单位m/s；
Vx0=0.12; %输入切向速度，单位m/s
Vx_ave=0.075; %输入平均切向速度，单位m/s
qu=1150*10^6; %输入材料屈服强度，单位Mpa，即：牛顿/平方米
P_ration=[0.3,0.22]; %泊松比，第一项为被加工工件的，第二项为磨粒的
E=[206*10^9,350*10^9]; %弹性模量（被加工工件、磨粒）
R=2*10^(-3); %磨粒半径
A_desnity=3850; %磨粒密度
HV1=680*10^3; %工件的维氏硬度（Gpa）；
KIC=130*10^4; %工件的断裂韧度
m=A_desnity*(4/3)*pi*R^3;  %单个磨粒的质量
V_A=(4/3)*pi*R^3; %计算单个磨粒的体积
Pt=654*10^6*9.8; %定义磨粒在切向方向收到的压力pa，即N/m^2
c3=0.95; %计算单个磨粒横向切割距离时用到的常数，试验测得
n=1.05; %切向方向上磨粒与工件接触面积的指数，用于计算Fx
def=zeros(numx,numy); %定义变形量矩阵
delta_h=zeros(numx,numy); %定义每个峰的去除高度矩阵
V_desnity=0.2624; %输入磨粒的体积浓度（每立方米）
T=180; %抛光时间（分钟）

%%计算法向最大变形深度

%弹性-塑性变形临界速度纵向Vy-ela-pla
Vy_ela_pla=(0.8*0.25^(5/2))^0.5*...
    ((qu^5/A_desnity)*((1-P_ration(1)^2)/E(1)+(1-P_ration(2)^2)/E(2))^4)^0.5*pi*pi;

%弹性-塑性变形临界变形深度
deform_ela_pla=0.25*pi^2*R*qu*((1-P_ration(1)^2)/E(1)+(1-P_ration(2)^2)/E(2))^2;
%弹性-塑性变形临界作用力
F_ela_pla=(2/3)*pi*qu*deform_ela_pla*R;
%塑性-脆性临界变形深度
deform_pla_bri=0.15*(E(1)/HV1)*(KIC/HV1)^2;
%塑性-脆性临界冲击速度
Vy_pla_bri=(Vy_ela_pla^2+qu*(3*deform_pla_bri+deform_ela_pla)*(deform_pla_bri-deform_ela_pla)/(4*R^2*A_desnity))^0.5;

%%计算变形深度（非永久/永久变形）deform/deform_per
if Vy0<=Vy_ela_pla  %弹性阶段

    deform=(0.8*pi*A_desnity*((1-P_ration(1)^2)/E(1)+(1-P_ration(2)^2)/E(2)))^0.4*R*Vy0^0.8;
    deform_per=0; %弹性阶段被加工工件表面不发生永久变形

end

if Vy0>Vy_ela_pla && Vy0<=Vy_pla_bri %塑性阶段
 
    %单个磨粒撞击工件表面产生撞击坑的最大深度
    deform=(1/3)*deform_ela_pla+(1/6)*...
        (16*deform_ela_pla^2+48*R^2*A_desnity*(Vy0^2-Vy_ela_pla^2)*qu^(-1))^0.5;  
    deform_per=deform-deform_ela_pla; %塑性阶段被加工工件表面产生的永久变形深度

end

%%计算单个磨粒切向方向的切割长度
ymax=deform; %宽度方向上被切割截面的最大深度
Dmax=2*sqrt(R^2-(R-deform)^2); %撞击坑最大直径，sqrt(R*deform)为撞击坑的半径
Ax=0.5*ymax*Dmax;
%%c3和n均为实验数据，需要手动输入

%L为单个磨粒的切割长度
L=0.5*c3*Vx0^2/(R^(1.25*n)*ymax^(0.75*n)*Pt); %c1,c2或c3均为实验测得的数据

l_desnity=(3*V_desnity/(4*pi*R^3))^(1/3); %磨粒的线密度
ls=1/l_desnity; %两磨粒的平均间距
Nx=(Lx+L)/ls; %取样区域内某一长度方向上的磨粒数
Ny=(Ly+2*R)/ls; %取样取样内某一宽度方向上的磨粒数
%N=Lx*Ly;
%N=(Lx*Ly)/(ls^2); %单位面积A=Lx*Ly内的磨粒数
N=Nx*Ny;

k = 1;
t_one=k*(ls/Vx_ave)/(60*N); %该表面每次被切割的时间间隔(分钟)
num=T/t_one; %表面被切割的次数

%被第num个磨粒作用后表面形貌的变化
for numm=1:1:num/2

    %%计算单个磨粒作用后的表面形貌变化
    i=randi([-floor(Dmax/(2*delta_dg))+2 numx+floor(Dmax/(2*delta_dg))]);
    j=randi([-floor(ls/(2*delta_dg))+2 numy+floor(ls/(2*delta_dg))]); %随机获得某个磨粒切割轨迹的中点

    nx=L/delta_dg; %单个磨粒在横向上切去的表面突起（峰）的个数
    z_mean=(sum(sum(z))/(numx*numy)); %计算表面峰高的平均值

    if deform_per==0
        def=zeros(numx,numy);   %若永久变形深度0，不产生形貌变化
    end

    %%

    if deform_per>0
        nxx=floor(L/delta_dg)+1; %横向上被切割过的峰的位置数量

        coutCutNum = 0;  %被切掉的峰的个数
        totolHigh = 0;  %背切的峰的高度和
       
        for nx=0:-1:-floor(L/(2*delta_dg))+1 %切向切割截面的左半部分

            if j+nx>0 && j+nx<=numy %判断在长度方向上是否超出取样区域

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                for ny=0:-1:round(-floor(D/(2*delta_dg))+2)

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %如果未超出范围，就记入到未切割的峰中
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end
            end
        end

        for nx=1:1:floor(L/(2*delta_dg))

            if j+nx>0 && j+nx<=numy

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                %对宽度方向上被切割到的点进行for循环迭代，以计算变形量
                for ny=0:-1:-D/(2*delta_dg)+1

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end
            end
        end

        mean_z = totolHigh/coutCutNum;
        
    %%
        for nx=0:-1:-floor(L/(2*delta_dg))+1 %切向切割截面的左半部分

            if j+nx>0 && j+nx<=numy %判断在长度方向上是否超出取样区域

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                %对宽度方向上被切割到的点进行for循环迭代，以计算变形量
                for ny=0:-1:round(-floor(D/(2*delta_dg))+2)

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end
            end
        end

        for nx=1:1:floor(L/(2*delta_dg))
            nxx=floor(L/delta_dg)+1; %横向上被切割过的峰的位置数量

            if j+nx>0 && j+nx<=numy

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                %对宽度方向上被切割到的点进行for循环迭代，以计算变形量
                for ny=0:-1:-D/(2*delta_dg)+1

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end
            end
        end
    end

    %%


    z=z-delta_h;
    def=zeros(numx,numy);
    delta_h=zeros(numx,numy);

end

z=z';

for numm=1:1:num/2

    %%计算单个磨粒作用后的表面形貌变化
    i=randi([-floor(Dmax/(2*delta_dg))+2 numx+floor(Dmax/(2*delta_dg))]);
    j=randi([-floor(ls/(2*delta_dg))+2 numy+floor(ls/(2*delta_dg))]); %随机获得某个磨粒切割轨迹的中点

    nx=L/delta_dg; %单个磨粒在横向上切去的表面突起（峰）的个数
    z_mean=(sum(sum(z))/(numx*numy)); %计算表面峰高的平均值

    if deform_per==0
        def=zeros(numx,numy);   %若永久变形深度0，不产生形貌变化
    end

    %%

    if deform_per>0
        nxx=floor(L/delta_dg)+1; %横向上被切割过的峰的位置数量

        coutCutNum = 0;  %被切掉的峰的个数
        totolHigh = 0;  %背切的峰的高度和
       
        for nx=0:-1:-floor(L/(2*delta_dg))+1 %切向切割截面的左半部分

            if j+nx>0 && j+nx<=numy %判断在长度方向上是否超出取样区域

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                for ny=0:-1:round(-floor(D/(2*delta_dg))+2)

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %如果未超出范围，就记入到未切割的峰中
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end
            end
        end

        for nx=1:1:floor(L/(2*delta_dg))

            if j+nx>0 && j+nx<=numy

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                %对宽度方向上被切割到的点进行for循环迭代，以计算变形量
                for ny=0:-1:-D/(2*delta_dg)+1

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        coutCutNum = coutCutNum + 1;
                        totolHigh = totolHigh + z(i+ny, j+nx);
                    end

                end
            end
        end

        mean_z = totolHigh/coutCutNum;
        
    %%
        for nx=0:-1:-floor(L/(2*delta_dg))+1 %切向切割截面的左半部分

            if j+nx>0 && j+nx<=numy %判断在长度方向上是否超出取样区域

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                %对宽度方向上被切割到的点进行for循环迭代，以计算变形量
                for ny=0:-1:round(-floor(D/(2*delta_dg))+2)

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end
            end
        end

        for nx=1:1:floor(L/(2*delta_dg))
            nxx=floor(L/delta_dg)+1; %横向上被切割过的峰的位置数量

            if j+nx>0 && j+nx<=numy

                def_central=((L-2*abs(nx)*delta_dg)/L)*deform_per; %宽度方向上中点处变形深度
                D=2*sqrt(R^2-(R-def_central)^2); %最大切割宽度

                %对宽度方向上被切割到的点进行for循环迭代，以计算变形量
                for ny=0:-1:-D/(2*delta_dg)+1

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end

                for ny=1:1:floor(D/(2*delta_dg))

                    if i+ny>0 && i+ny<=numx  %判读在宽度方向上是否超出取样范围
                        %计算点（i+nx，j+ny）处的变形量
                        def(i+ny,j+nx)=def_central-R+sqrt(R^2-(nx*delta_dg)^2);
                        delta_h(i+ny,j+nx)=def(i+ny,j+nx)*10^6-mean_z+z(i+ny,j+nx);
                        if delta_h(i+ny,j+nx)<0
                            delta_h(i+ny,j+nx)=0;
                        end
                    end

                end
            end
        end
    end

    %%


    z=z-delta_h;
    def=zeros(numx,numy);
    delta_h=zeros(numx,numy);

end



%%计算工件抛光后的表面粗糙度
z_ave=sum(sum(z))/(numx*numy);
z_norm=z-z_ave*ones(numx,numy);
zabs_polish=abs(z_norm); 
Ra_polish=double(sum(sum(zabs_polish))/(numx*numy));




%被加工之后的表面形貌figure(2)
figure(2)
surf(xx1,yy1,z_norm,'EdgeColor','interp');
%surf(xx1,yy1,z,'EdgeColor','interp');
%zlim([-2 2]);
view(-60,80);
%title('模型预测形貌图');
colorbar;